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In this paper we present a mathematical algorithm for constructing 
a smoothing cubic spline with periodic end conditions and a prede- 
termined 'closeness of fit' to a given set of points in the plane. In 
addition to providing a mathematical tool for smoothing raw data in 
which the underlying function is known to be periodic, this algorithm 
has special significance in computer graphics, because the use of 
smoothing functions with periodic end conditions is essential for 
producing visually acceptable, smooth, closed curves. Sample plots 
are included to illustrate the power and flexibility of this algorithm. 

I. INTRODUCTION 

Although "natural" splines are used extensively and are quite ap- 
propriate for smoothing many types of data, they often produce less 
than satisfactory results when used to smooth data points that belong 
to a periodic function. The inappropriateness of using "natural" splines 
to approximate periodic data is especially evident in graphics appli- 
cations. In particular, when the data points represent a closed curve, 
smoothing (parametricaJly) with "natural" splines will lead to unac- 
ceptable results because the "natural" end conditions will cause the 
curve either to close up with a noticeable discontinuity, or to not close 
up at all (see Fig. 1). 

Existing methods for constructing smoothing splines with a prede- 
termined closeness of fit all lead to splines with "natural" end condi- 
tions. 1,2 A method developed by Spath 3 produces a smoothing spline 
with periodic end conditions, but the closeness of fit cannot be deter- 
mined in advance. In this paper we will describe a method for con- 
structing a smoothing cubic spline that has periodic end conditions 
and that also satisfies a predetermined closeness of fit to a given set of 
data points. 

This algorithm has potentially wide applicability, especially in the 
realm of interactive graphics. It makes possible the computer genera- 
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Fig. 1 — Comparison of periodic versus natural spline smoothing, (a) Periodic cubic 
spline smoothing with uniform weights, (b) Natural cubic spline smoothing with small 
weights at end points, (c) Natural cubic spline smoothing with uniform weights. 

tion of free-form, smooth, closed curves by merely specifying the 
approximate locations of as few as three distinct points. The shape of 
the curve can be controlled easily by moving one or more points, or by 
adjusting the weighting factors associated with some or all of the 
points. 

An efficient program based on this new algorithm has been written 
and tested. Sample plots illustrating this method are included. 

II. TERMINOLOGY 

Let Pk = (xk, yk), k = 1, n, be n points in the plane. A "cubic spline" 
on [x\, x n ] with knots at X\, • • • , x n , is a function /"that coincides with 
a third-order polynomial fk on each sub-interval [xk, Xk+i], k = 1, n — 
1, and such that / is continuous and has continuous first and second 
derivatives over the entire interval [xi, x n ]. 

In other words, / is a cubic spline on [xi, x n ] if, for each k = 1, n — 
1 there exist real numbers a«, bk, C/ e , dk (the "spline coefficients" of/ - ) 
such that for every x in [xk, Xk+i], 
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f( x ) = f k ( x ) = a k + b k {x - x k ) + c k (x - x k ) 2 + d k (x - x*) 8 . (1) 

Furthermore, the continuity off, f, and f" on [xi, x„] implies that at 
each interior knot Xk, k = 2, n — 1, 

fk-i(x k ) = fkix*), (2) 

f'k-x(xk) = ftixk), (3) 

fk-i(Xk) - fk (Xk). (4) 

The cubic spline f is said to be "periodic" if it satisfies the following 
additional conditions (known as "periodic end conditions"): 

f(Xn)=f(Xl), (5) 

f'(Xn)=f'(X 1 ), (6) 

f"(x n )=f"( Xl ). (7) 

A "natural" cubic spline differs from a "periodic" cubic spline in that 
it satisfies the so-called "natural end conditions": f"{x x ) = f"(x„) = 0. 

III. FORMAL STATEMENT OF THE PROBLEM 

Let Pa- = (x/,, yk),k = 1, n, be n points in the plane, with x x < x 2 < 
• •• < x n . Let Wk, k = 1, n, be positive real numbers ("weighting 
factors") associated with P k , k = 1, n, respectively. (Assume ;y„ = y it 
w n = Wi.) Given an arbitrary constant, M > 0, the problem is to 
determine the set of 4(n — 1) coefficients of the periodic cubic spline 
/"on [x\ , x n ] with knots at Xi , • • • , x n , such that 

(i) f has minimal "total curvature" G(f) = /*; f"{xfdx, and 
{ii) f satisfies the following weighted, distance-squared constraint, 
or "closeness of fit," with respect to the given points: 



mn= i 



*-] 



f(xk) - yk 



Wk 



<M. 



Note that the weighting factors give inverse importance to the 
points. Note also that M controls the degree of smoothness, so that 
increasing the value of M for a fixed set of weighting factors will lead 
to a smoother, or flatter, spline. Conversely, choosing a sufficiently 
small value of M will lead to a spline that closely approximates an 
interpolating spline. In general, an appropriate choice for the value of 
M will depend on the values chosen for the weighting factors. If, for 
example, the weighting factors are chosen so that W/, is the standard 
deviation at x*, then a suitable choice for M would be some value 
between the confidence limits, n — V2n and n + V2n. 
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IV. GENERAL APPROACH TO THE SOLUTION 

To minimize G{ f) subject to the constraint H{ f) < M, we introduce 
an auxiliary variable z and a non-negative Lagrange multiplier p and 
minimize 

F(f) = G(f)+p[H(f) + z 2 -M]. 

This is the approach used by Reinsch in Ref. 2 to solve the analogous 
problem for "natural" cubic splines. However, Reinsch minimizes F 
over the class of all continuous functions with continuous first and 
second derivatives on [xi, x„], which leads to a "natural" cubic spline 
as the solution. We will restrict the class of admissible functions in the 
minimization of Fto periodic cubic spline splines on [xi, x„] with knots 
at Xi, • • • , x„ , obtaining a direct solution to our problem. 

V. LINEAR SYSTEM RESULTING FROM CONTINUITY AND PERIODICITY 
CONDITIONS 

Let / be an arbitrary periodic cubic spline on [xi, x„], with spline 
coefficients a*, bk, Ck, d*» k ™ 1, n — 1. Then, for x in [x/,, Xk+\], f(x) is 
given by (1) above, and f'(x) and f"(x) are given below: 

f'(x) = fUx) = b k + 2c k (x - x k ) + 3d k (x - x k )\ (8) 

f" (x) = fS (x) = 2c k + 6d k (x - x k ). (9) 

Expressing f, f, and f" explicitly in eqs. (1), (8), and (9), respectively, 
allows us to derive linear relationships among the spline coefficients. 
From the continuity off at the interior knots of (3) and the periodicity 
of f on [X] , x n ] in (6), it follows that: 

2c k h k = bk + i - b k - 3dM, for k = 1, n - 1, (10) 

where hk = x/,+i — Xk for k ■ 1, n — 1, and b n denotes b\. 

From the continuity of /at the interior knots of (2) and the period- 
icity of /on [xi , x„] in (5), the first-order coefficients may be expressed 
in terms of the constant, second-order, and third-order coefficients: 

bk = (a*+i - a k )/hk - c k h k - d k h\, for k = 1, n - 1, (11) 

where a„ denotes a.\ . 

From the continuity of /" at the interior knots of (4) and the 
periodicity of/" on [x i} x n ] in (7), the third-order coefficients may be 
expressed as a function of the second-order coefficients: 

d k = (c k+ i - c k )/3hk, for k = 1, n - 1, (12) 

where c„ denotes C\ . 

Using (11) and (12) to eliminate the bk 's and dk's from (10) leads to 
a system of linear equations in the a.k's and c*'s, given in matrix 
notation as follows: 
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Sc = 3Qa, (13) 

where S and Q are symmetric, cyclic-tridiagonal matrices of order n 
— 1, and c, a are the column vectors (ci, • • • , c n -i) , (ai, ••• , a n -\) , 
respectively. 

The non-zero entries of S and Q are expressed in terms of the 
distances between successive knots: 

S(A>, k) = 2(h k -i + h k ) for k = 1, n - 1 
S(k, k + l) = S{k + 1, k) = h k for k - 1, n - 2 
8(1, n - 1) = S(n - 1, 1) = ft„-i 

Q(&, &) = -l//i A -i - \/h k for ft = 1, n - 1 

Q(£, ^ + 1) = Q(A: + 1, k) = 1/hk for k = 1, n - 2 

Q(l, » - 1) = Q(n - 1, 1) = l/hn-i, 

where /i denotes h n -\ . By Gershgorin's Theorem 4 it can be shown that 
S is positive definite (and therefore non-singular), while Q is positive 
semi-definite and singular with rank n — 2. 

VI. LINEAR SYSTEM RESULTING FROM MINIMIZING F WITH RESPECT 
TO THE CONSTANT COEFFICIENTS 

Note that (13) is a system of n - 1 linear equations in 2(n - 1) 
unknowns: the constant coefficients and the second-order coefficients. 
We shall derive a second system of n — 1 linear equations in these 
unknowns. We proceed by first showing that H and G can be expressed 
as functions of the constant coefficients only. 

From the spline representation of / in (1), it follows immediately 
that 






f{x k ) - y k 



2 n / V 2 

= y ( ak-y >< 



Wk 

is a function of a x , • • • , a„. And since the periodicity of /"implies 

a„ =f(x n ) = f(xi) = a u 

His a function of the constant spline coefficients a it • • • , a n -\. 

From the explicit representation of f" in (9), the "total curvature" 
G can be expressed as 

f"{x) 2 dx=l f"(x) 2 dx 

= X [2c k + 6d k (x - x k ) 2 dx. 

*-» Jr. 
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Evaluating the integrals over each sub-interval directly and eliminating 
the dk's with (12) lead to the following: 
«-i 4 
G = X - h k {cl + c k c k+ i + cl+i). 

k=l o 

Rewriting in matrix notation and applying (13), we have: 

G = (%)c r Sc = (%)(3S- 1 Qa)S(3S- 1 Qa) = Ga^S'Qa. 

Thus, G also can be expressed in terms of the constant spline coeffi- 
cients. Note that from this representation of G, 

fi=12Q,S- , Qa 
dak 

for k = 1, n — 1, where Q* is the £th row of Q. 

Since G and H are functions of a\ , • • • , a n -\ , then F is a function of 
the independent variables Oi, • • • , a n -\,p, and z. In order for F to be 
minimized, the partial derivative of F with respect to each of its 
independent variables must vanish. Thus, for each k = 1, n — 1, 
differentiating F with respect to a* yields: 

f^f + « 12Q4 s-Qa + 2p (f^ 

dan oak dak \ w% 



Rewriting this set of n — 1 linear equations in matrix notation and 
using (13) to replace S _1 Qa with c/3 lead to: 

4Qc + 2pW- 2 (a - y) = 0, (14) 

where y is the column vector (yi, • • • , y n -\) T . Combining (13) and (14), 
we have the following linear system in c: 

( P S + 6QW 2 Q)c = 3pQy. (15) 

The matrix A p = {pS + 6QW 2 Q) is symmetric, five-banded with three 
non-zero entries in the upper right and lower left corners. It can be 
shown that, for all positive values of p, A p is positive definite. Thus, 
for each/? > 0, (15) has a unique solution in c. 

Note that for each non-zero value of p, a = y — (2//?)W 2 Qc from 
(14). Note also that from (12) d is uniquely determined by c, and from 
(11) b is uniquely determined by a, c, and d. Thus, to each positive 
value of p corresponds a unique periodic cubic spline on [xi, x n ], whose 
coefficients are given in the vectors a, b, c, d. 



VII. CONSEQUENCES OF MINIMIZING F WITH RESPECT TO p AND z 

Minimizing F = G + p(H + z 2 — M) with respect to the Lagrange 
multiplier p leads to: 
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BF 

— = H+z 2 -M = 0, 

dp 

which merely states that the distance constraint on H (expressed as 
an equality in terms of the auxiliary variable z) must be satisfied when 
the minimal value of F is attained. 
On the other hand, minimizing F with respect to z yields: 

dz 

which implies that at least one of the two variables, p or z, must be 
equal to when the minimal value of F is attained. 

Note that if p = when F is minimized, then F = G. Since G = 
(2/3)c r Sc and S is positive definite, then G is minimized when c = 0. 
This in turn implies d = 0, so that the second- and third-order 
coefficients vanish, resulting in a piecewise linear rninimizing spline. 
The properties of being piecewise linear and having a continuous first 
derivative together imply that the minimizing spline is a straight line. 
Furthermore, periodicity of the spline implies that the straight line is 
in fact horizontal. 

On the other hand, if p > when F is minimized, then z = 0, so that 
H = M. Since a = y - (2/p)W 2 Qc from (14), and c = 3pA^Qy from 
(15), if can be expressed as a function of p. Thus, if minimization of F 
occurs for a positive value of p, it remains to determine the value of p 
for which H(p) = M. 

VIII. PROPERTIES OF H AS A FUNCTION OF p 

The following facts can be established: for all positive values of p, 
H(p) is a continuous, convex function of p with negative slope. Fur- 
thermore, as p approaches zero from the right, H(p) becomes arbi- 
trarily large. 

IX. ALGORITHM FOR DETERMINING SPLINE COEFFICIENTS 

We can now state the following algorithm for determining the 
minimizing spline. Compute the equation of the horizontal straight 
line with the least-squares fit to the given data points: 

\ k % wljl \*-i wi) 

Determine if this line satisfies the distance constraint on H(f). If it 
does, we are done. If it does not satisfy the distance constraint, start 
with some positive value of p and search for the value of p for which 
H(p) = M, using a combination of Newton's method when moving to 
the right and a binary search when moving to the left (or any applicable 
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search). Insert this value of p in (15), solve the linear system for c, and 
compute the related values of a, b, and d using (14), (11), and (12). 
The periodic cubic spline associated with this value of p will be our 
solution. 

X. PRACTICAL CONSIDERATIONS IN SOLVING THE LINEAR SYSTEM 

Since the matrix A p is positive definite for each p > 0, it can be 
decomposed into the product of a lower triangular matrix R and its 
transpose, using the square-root (Cholesky's) method. 5 The linear 
system A p c = 3pQy can then be solved efficiently in two steps, by 
applying forward substitution to the lower triangular system Rv = 
3/?Qy, followed by backward substitution to the upper triangular 
system R T c = v. 

Furthermore, since A p is symmetric and five-banded with three non- 
zero entries in two corners, its decomposition R will consist of three 
non-zero bands (the main diagonal and the two diagonals below it) 
and two non-zero rows along the bottom, so that R can be stored in 
fewer than bn locations. (The entries of A p need not be stored; they 
may be computed as needed.) 

The sparseness of the matrix R and its structure described above 
lead not only to its compact storage, but also to the linear time solution 
of the upper and lower triangular systems, and hence to the linear time 
solution of the system A p c = 3pQy, for each non-zero value of p. 

It should be pointed out that, unless the number of data points to be 
smoothed is rather limited (approximately 30 or fewer), the straight- 
forward application of Cholesky's method to decompose A p will en- 
counter underflow problems. This is due to the fact that, as the 
dimension of A p increases, entries with exponentially decreasing mag- 
nitudes will appear in its decomposition R. This difficulty can be 
circumvented by truncating sufficiently small values to zero, while still 
retaining single-precision accuracy in the solution of the triangular 
systems. (Truncation has the additional advantage of significantly 
reducing the number of arithmetic operations required in computing 
the entries of R when the number of data points is large.) 

The program implementing this algorithm has been tested on an 
IBM-370 computer using single-precision arithmetic, and has success- 
fully smoothed up to 250 points before encountering detectable round- 
off errors. On the average, the Newton and binary search converged 
after six to eight iterations, independent of the number of data points. 

XI. SAMPLE PLOTS 

The data points in Figs. 1 and 2 were generated by adding random 
noise to an ellipse. The dotted curves represent cubic spline interpo- 
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Fig. 2— Periodic cubic spline smoothing for a varying number of data points with 
uniform weights. 

lation of the data (using periodic splines), while the solid curves 
represent cubic spline smoothing of the same data. In each case a 
parameter was introduced so that the curve could be represented as 
two separate single- valued functions of the parameter. Then smoothing 
was performed twice, once with the x values as a function of the 
parameter and then again with the y values as a function of the 
parameter. The smoothed x and smoothed v values were then plotted 
against the parameter to produce the closed curve. Figure 2 illustrates 
the algorithm with a varying number of data points, from only three 
distinct points to 250 points. In each case, uniform weights were used. 
A "tight" fit was chosen in the example with three points to show how 
the method can be used to simulate periodic interpolation of the 
points. 
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